# %% [markdown]
""" 
# 案例数据背景介绍
"""
# %% [markdown]
"""
1.数据来源：

    Dongjin Cho，韩国蔚山国立科技学院城市与环境工程学院

    Cheolhee Yoo，韩国蔚山国立科技学院城市与环境工程学院
 """
# %% [markdown]
# 2.摘要：
#
#     本案例所使用的数据集包含了韩国首尔2013至2017年夏天的14个天气预报（NWP）的气象预报数据、2个现场观测数据和5个地理辅助变量。
#
#
# %% [markdown]
# 3.数据集信息：
#
#     这些数据是为了修正韩国气象局在首尔运行的LDAPS(Local Data Assimilation and Prediction System)模式下次日最高和最低气温预报的偏差。该数据包括2013年至2017年的夏季数据。输入数据主要由LDAPS模型的次日预报数据、当前的现场最高和最低气温以及地理辅助变量组成。此数据有两个输出数据（即第二天的最高和最低气温）；2015年至2017年期间进行了后期验证。
#     关于LDAPS模式: http://qikan.camscma.cn/jamsweb/cn/article/doi/10.11898/1001-7313.20200103?viewType=HTML
#
# %% [markdown]
"""
4.数据属性介绍：
 
       1. station - 气象站编号: 1 - 25 
       
       2. Date - 记录日期: yyyy-mm-dd ('2013-06-30' to '2017-08-30') 
       
       3. Present_Tmax - 记录日期的0-21时最高气温(Â°C): 20 - 37.6 
       
       4. Present_Tmin - 记录日期的0-21时最低气温(Â°C): 11.3 - 29.9 
       
       5. LDAPS_RHmin - LDAPSM模式预测的次日最小相对湿度 (%): 19.8 - 98.5 
       
       6. LDAPS_RHmax - LDAPS模式预测的次日最大相对湿度 (%): 58.9 - 100 
       
       7. LDAPS_Tmax_lapse - LDAPS模式预测的次日最高气温递减率 (Â°C): 17.6 - 38.5 
       
       8. LDAPS_Tmin_lapse - LDAPS模式预测的次日最高气温递减率 (Â°C): 14.3 - 29.6 
       
       9. LDAPS_WS - LDAPS 模式预测的次日平均风速 (m/s): 2.9 to 21.9 
       
       10. LDAPS_LH - LDAPS模式预测的次日平均潜热通量 (W/m2): -13.6 - 213.4 
       
       11. LDAPS_CC1 - LDAPS模式预测的次日第一个6小时平均云量 (0-5 h) (%): 0 - 0.97 
       
       12. LDAPS_CC2 - LDAPS模式预测的次日第二个6小时平均云量 (0-5 h) (6-11 h) (%): 0 to 0.97 
       
       13. LDAPS_CC3 - LDAPS模式预测的次日第三个6小时平均云量 (0-5 h) (12-17 h) (%): 0 to 0.98 
       
       14. LDAPS_CC4 - LDAPS模式预测的次日第四个6小时平均云量 (0-5 h) (18-23 h) (%): 0 to 0.97 
       
       15. LDAPS_PPT1 - LDAPS模式预测的次日第一个6小时平均降水量 (0-5 h) (%): 0 to 23.7 
       
       16. LDAPS_PPT2 - LDAPS模式预测的次日第二个6小时平均降水量  (6-11 h) (%): 0 to 21.6 
       
       17. LDAPS_PPT3 - LDAPS模式预测的次日第三个6小时平均降水量 (12-17 h) (%): 0 to 15.8 
       
       18. LDAPS_PPT4 - LDAPS模式预测的次日第四个6小时平均降水量  (18-23 h) (%): 0 to 16.7 
       
       19. lat - 维度 (Â°): 37.456 - 37.645 
       
       20. lon - 经度 (Â°): 126.826 - 127.135 
       
       21. DEM - 高程 (m): 12.4 - 212.3 
       
       22. Slope - 坡度 (Â°): 0.1 - 5.2 
       
       23. Solar radiation - 太阳辐射量 (wh/m2): 4329.5 to 5992.9 
       
       24. Next_Tmax - 次日最高气温 (Â°C): 17.4 to 38.9 
       
       25. Next_Tmin - 次日最低气温 (Â°C): 11.3 to 29.8 """


# %% [markdown]
# **5. sklearn中常见的几个线性回归类库：**
#
# ![image.png](attachment:image.png)
# %% [markdown]
#
#
#
# %% [markdown]
# # 读取数据

# %%
import numpy as np
import pandas as pd
import pandas_profiling as pp
import matplotlib.pyplot as plt
import seaborn as sns

# 全部行都能输出
from IPython.core.interactiveshell import InteractiveShell

InteractiveShell.ast_node_interactivity = "all"


# %%
Data = pd.read_csv(r"Dataset.csv")
Data.head()

# %% [markdown]
# # 数据探索
# %% [markdown]
# ## 基本统计信息

# %%
Data.describe().T


# %%
Data.info()

# %% [markdown]
# ## 查看空值（5分）

# %%
Data.isna().sum()


# %%
# 取出所有有空值的记录  （5分）
Data[Data[["Next_Tmax", "Next_Tmin"]].isna().sum(axis=1) != 0].shape
# %% [markdown]
# ## 查看每列数据的分布状态

# %%
# 绘制每列数据的直方图（分布频次）
plot = Data.drop("Date", axis=1)  # 删除'Date'列

for I in plot.columns:
    sns.distplot(plot[I])
    plt.show()

# %% [markdown]
# ## 生成数据概况报告

# %%
# 生成整体数据概要
data_report = pp.ProfileReport(Data)
# data_report

# 也可以把生成的报告导出保存
data_report.to_file("report.html")

# %% [markdown]
# # 数据预处理
# %% [markdown]
# ## 提取数据

# %%
# 由于本数据集的特殊性（有两个样本标签：次日最高气温和最低气温）
# ，本案例中我们只进行次日最高温预测--因此这里只选出与最高温预测相关的变量，作为建模要用的数据

Max = Data[
    [
        "Present_Tmax",
        "LDAPS_RHmax",
        "LDAPS_Tmax_lapse",
        "LDAPS_WS",
        "LDAPS_LH",
        "LDAPS_CC1",
        "LDAPS_CC2",
        "LDAPS_CC3",
        "LDAPS_CC4",
        "LDAPS_PPT1",
        "LDAPS_PPT4",
        "lat",
        "lon",
        "DEM",
        "Slope",
        "Solar radiation",
        "Next_Tmax",
    ]
]


# %%
Max.head()


# %%
Max.info()


# %%
Max.describe().T

# %% [markdown]
# ## 填充缺失值（10分）

# %%
# 用每列的均值填充其缺失值
for c in Max.columns:
    m = Max[c].mean()
    print(m)
    Max[c] = Max[c].fillna(m)

# %%
Max.isnull().sum()

# %% [markdown]
# ## 标准化处理连续型特征

# %%
from sklearn.preprocessing import StandardScaler

# 标准化处理连续型特征 （10分）

# %% [markdown]
# ## 划分特征列和标签列
Scaled_Data = pd.DataFrame(StandardScaler().fit_transform(Max), columns=Max.columns)
Scaled_Data
# %%
X = Scaled_Data.drop(["Next_Tmax"], axis=1)
X.head()


# %%
Y = Scaled_Data["Next_Tmax"]
Y

# %% [markdown]
# # 创建线性回归模型并评估
# %% [markdown]
# ## 划分数据集（5分）

# %%
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, Y, test_size=0.33, random_state=43
)
# %% [markdown]
# ## 建模（5分）
from sklearn.linear_model import LinearRegression


# %%
LR = LinearRegression().fit(X_train, y_train)

# %% [markdown]
# ## 评估（5分）
LR.score(X_train, y_train), LR.score(X_test, y_test)

# %% [markdown]
# ## 计算MSE与RMSE（10分）
from sklearn import metrics

# %%
print(
    """ 
    mse:{}
    rmse:{}
    """.format(
        metrics.mean_squared_error(y_test, LR.predict(X_test)),
        metrics.mean_squared_error(y_test, LR.predict(X_test)) ** (1 / 2),
    )
)


# %% [markdown]
# # 模型诊断
# %% [markdown]
# ## 绘制线性拟合图（10分）
# %% [markdown]
# '''
# 从上可以看到RMSE(越小越好）的值0.5，我们通过可视化来看一下训练后的预测和真实值之间的差异：
# '''

# %%
# 画个折线图来看看模型拟合的好坏程度
y_pred = LR.predict(X_test)
plt.figure(figsize=(20, 10))
plt.plot(range(len(y_test)), y_test, "r", label="Y-TEST")
plt.plot(range(len(y_test)), y_pred, "b", label="Y-PREDICTED")
plt.legend()

# %% [markdown]
# '''
# 可以看到许多地方，红线会明显超出蓝线，说明我们的模型拟合度不够，再换个散点图来更直观地看一下：
# '''
# %%
# 使用散点图观测（10分）

plt.figure(figsize=(12, 12))
plt.scatter(y_test, y_pred)
l = y_test
plt.plot(y_test, l, c="black", linestyle="--")
plt.xlabel("Y-TEST")
plt.ylabel("Y-PREDICTED")

# %% [markdown]
# '''
# 图中我们可以看到，如果完全拟合，散点应该和直线相重合，这里发现，y_test有一些的异常值，
# 而线性回归模型的一大缺点就是对异常值很敏感，会极大影响模型的准确性，因此，下一步，我们就根据这一点，对模型进行优化
# '''
# %%
# 绘制残差图查看模型的拟合情况（10分）
from sklearn.model_selection import cross_val_predict, cross_val_score

# %% [markdown]
# ## 绘制残差图（10分）
#%%
plt.figure(figsize=(15, 5))
plt.scatter(Y, Y - LR.predict(X), s=1)
plt.axhline(0)


# %% [markdown]
# ## 绘制相关性热图（10分）
#
# **绘制各预测变量间的相关性热图，查看是否存在相关性很高的变量（共线性或近似共线性）**

# %%
plt.figure(figsize=(15, 17))
corrh = X_train.corr()
plt.imshow(corrh, cmap=plt.cm.hot)
plt.xticks([i for i in range(len(X_train.columns))], X_train.columns, rotation=90)
plt.yticks([i for i in range(len(X_train.columns))], X_train.columns)
for i in range(corrh.shape[0]):
    for j in range(corrh.shape[1]):
        plt.text(
            j,
            i,
            corrh.iloc[i, j].round(2),
            ha="center",
            va="center",
            color="w",
        )

plt.colorbar()
# %% [markdown]
# # 模型调优
# %% [markdown]
# ## 删除异常值（10分）

# %%
#####================寻找异常值====================######
Scaled_Data["Next_Tmax"].plot.box()


# %%
# 删除掉Next_Tmax中小于下边缘值的记录（10分）

Qs = Scaled_Data["Next_Tmax"].quantile(q=[0.25, 0.75])
low_whisker = Qs[0.25] - 1.5 * (Qs[0.75] - Qs[0.25])
low_whisker

# %%
# Y中小于箱型图下边缘的值
Scaled_Data[Scaled_Data["Next_Tmax"] < low_whisker]


# %%
# 要删除的样本的Index
drop_index = Scaled_Data[Scaled_Data["Next_Tmax"] < low_whisker].index.values
drop_index.shape


# %%
# 删除Y中小于箱型图下边缘的值
X = X.drop(drop_index)
Y = Y.drop(drop_index)


# %%
# 删除后的X,Y -------少了两行
X
Y


# %%
# 重新划分数据集
from sklearn import metrics

X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.3, random_state=1)
# 对训练集进行训练
lr = LinearRegression()
lr.fit(X_train, Y_train)
# 对测试集进行预测
Y_pred = lr.predict(X_test)
MSE = metrics.mean_squared_error(Y_test, Y_pred)
RMSE = np.sqrt(metrics.mean_squared_error(Y_test, Y_pred))
print("MSE:", MSE)
print("RMSE:", RMSE)


# %%
# 上一次评估结果
# MSE: 0.25266554606897074
# RMSE: 0.502658478560713

""" 
mse:0.2355301777756762
rmse:0.48531451428499045 """
# %% [markdown]
# '''
# 可以看到，该模型的RMSE由原来的0.5降低到了0.48，小幅度的提升了模型的准确率;去除异常值，提高准确度。
# '''
# %% [markdown]
# ## 特征筛选

# %%
X


# %%
# 特征筛选，选出和样本标签之间存在显著相关性（p<0.05)的特征，从而过滤掉对模型预测贡献不大的特征
from sklearn.feature_selection import f_regression

F, p = f_regression(X, Y)  # 一个特征对应一个F值和p值
k = F.shape[0] - (p > 0.05).sum()  # 去掉p值<=0.05的特征，k就是剩下的特征数量
k


# %%
# F检验筛选特征
from sklearn.feature_selection import SelectKBest

X_f = SelectKBest(f_regression, k=k).fit_transform(X, Y)
X_f.shape  # 过滤掉了两个特征
X_f


# %%
xtrain, xtest, ytrain, ytest = train_test_split(X_f, Y, test_size=0.3, random_state=1)


# %%
# # 实例化并训练模型
lr = LinearRegression().fit(xtrain, ytrain)


# %%
# 测试集上的评分（R²）
lr.score(xtrain, ytrain), lr.score(xtest, ytest)


# %%
ypred = lr.predict(xtest)
MSE = metrics.mean_squared_error(ytest, ypred)
RMSE = np.sqrt(metrics.mean_squared_error(ytest, ypred))
print("MSE:", MSE)
print("RMSE:", RMSE)

# %% [markdown]
# ## PCA降维
#
# 主成分分析的目的是从现有的特征中重建新的特征，新的特征剔除了原有特征中的冗余信息，因此更有区分度。
#
# 主成分分析的结果是得到新的特征，而不是简单地舍弃原来的特征列表中的一些特征。

# %%
from sklearn.decomposition import PCA

pca = PCA(n_components=0.97).fit(X, Y)
X_pca = pca.transform(X)


# %%
X_pca = pd.DataFrame(X_pca)
X_pca


# %%
Y


# %%
components = pca.components_  # 14个新特征的特征向量(每一行)
components = pd.DataFrame(components)
components


# %%
pca.explained_variance_ratio_  # 主成分的解释方差百分比


# %%
# 重新划分训练集和测试集
Xtrain, Xtest, Ytrain, Ytest = train_test_split(X_pca, Y, test_size=0.3, random_state=1)


# %%
# # 实例化并训练模型
lr = LinearRegression().fit(Xtrain, Ytrain)


# %%
# 模型在训练集和测试机上的评分（R²）
lr.score(Xtrain, Ytrain), lr.score(Xtrain, Ytrain)

# %% [markdown]
# # 随机森林回归 RandomForestRegression（20分）

# %%
# 1.使用集成算法-随机森林回归器
# 2.使用网格搜索寻找主要参数的最优取值组合
# 3.使用去除异常值后的训练集数据训练模型（15分）
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestRegressor

RF = RandomForestRegressor(criterion="mse")
P = ({"max_depth": range(1, 10), "n_estimators": range(1, 200, 10)},)
#%%
GS = GridSearchCV(
    cv=5,
    error_score=np.nan,
    estimator=RF,
    #  iid='deprecated',
    n_jobs=2,
    param_grid=P,
    pre_dispatch="2*n_jobs",
    refit=True,
    return_train_score=False,
    scoring=None,
    verbose=0,
)
GS.fit(X_train, Y_train)
# 查看最优参数组合
GS.best_params_


# %%
# 用上面得到的最优参数组合重新实例化模型
rfr = RandomForestRegressor(
    n_estimators=171,
    max_depth=9,
).fit(X_train, Y_train)


rfr.score(X_train, Y_train), rfr.score(
    X_test, Y_test
)  # 随机森林回归算法的score返回的是R²，并不是mse(尽管实例化时用的criterion='mse')


# %%
# # 使用交叉验证---随机森林回归器在测试集上的MSE平均得分（5分）
from sklearn.model_selection import cross_val_score

-cross_val_score(rfr, X_test, Y_test, cv=5, scoring="neg_mean_squared_error").mean()

# %% [markdown]
# '''
# 可以看出：无论是在训练集上还是测试集上，随机森林的的表现均优于LinearRegression,Lasson和Ridge：
# 训练集上的MSE由原来的0.25降到了0.15，R²由原来的0.75上升到了0.85。
# '''

# %%
